ACoP10 2019
Metrum Research Group
. # A tibble: 1 x 1
. Median_CFB
. <dbl>
. 1 -14.4
mrgsolveR package for simulation from ODE-based models
C++ inside model specification formatODEPACK / DLSODA (FORTRAN)
RID, \(\eta\), \(\varepsilon\)F, ALAG, SS etc, handled under the hoodR, C++, Rcpp, boost, RcppArmadilloR is it’s natural habitatmrgsolve started as QSP modeling toolR front end to deSolveC++ interface to DLSODARcpp - vectors, matrices, functions, environments, random numbersboost - numerical tools in C++C++ code (functions, data structures, classes)SBML to mrgsolve using R bindings to libSBMLGitHub site: https://github.com/metrumresearchgroup/mrgsolve
Issues and questions: https://github.com/metrumresearchgroup/mrgsolve/issues
mrgsolve website: https://mrgsolve.github.io
User Guide: https://mrgsolve.github.io/user_guide
Vignettes: https://mrgsolve.github.io/vignettes
The first / simplest workflow.
that come from?event in this exampleplot, mutate, as_data_frame etc …%>% %>% %>% Better.
%>% ...
----------------- source: pk1.cpp -----------------
project: /Users/kyleb/Rli...gsolve/models
shared object: pk1-so-4d02475ef035
time: start: 0 end: 192 delta: 0.2
add: <none>
compartments: EV CENT [2]
parameters: CL V KA [3]
captures: CP [1]
omega: 0x0
sigma: 0x0
solver: atol: 1e-08 rtol: 1e-08 maxsteps: 20k
------------------------------------------------------
Model parameters (N=3):
name value . name value
CL 1 | V 20
KA 1 | . .
CL, V2, Q etcTHETAs)WT, EGFR etc)
Model initial conditions (N=2):
name value . name value
CENT (2) 0 | EV (1) 0
Quiz:
. [1] "/Users/kyleb/Rlibs/mrgsolve/models"
.
.
. ---------------- source: effect.cpp ----------------
.
. project: /Users/kyleb/Rli...gsolve/models
. shared object: effect-so-4f5677552ab8
.
. time: start: 0 end: 36 delta: 0.1
. add: <none>
.
. compartments: GUT CENT PERIPH CE [4]
. parameters: VC KA K10 K12 K21 E0 EMAX EC50 KEO [9]
. captures: EFFECT CP [2]
. omega: 0x0
. sigma: 0x0
.
. solver: atol: 1e-08 rtol: 1e-08 maxsteps: 20k
. ------------------------------------------------------
. Error : the model file simple.cpp does not exist.
. [1] "Error : the model file simple.cpp does not exist.\n"
. attr(,"class")
. [1] "try-error"
. attr(,"condition")
. <simpleError: the model file simple.cpp does not exist.>
Summary:
mread to read, compile and load a modelWe use mread() to get a model loaded into R
Next, we’ll need
model %>% %>% simulate %>% take-a-look
Event object = quick / easy way to implement dose or other intervention
. Events:
. time amt cmt evid
. 1 0 100 1 1
time, evid, cmtev(...)evidevid 1, with rate==0)evid 1, with rate > 0)evid 2)
evid 3)evid 4)evid 8). time amt cmt evid
. 1 0 100 1 1
We will use a
model %>% intervention %>%
. Model: effect
. Dim: 362 x 8
. Time: 0 to 36
. ID: 1
. ID time GUT CENT PERIPH CE EFFECT CP
. 1: 1 0.0 0.00 0.000 0.0000 0.0000 157.0 0.000
. 2: 1 0.0 100.00 0.000 0.0000 0.0000 157.0 0.000
. 3: 1 0.1 91.21 8.443 0.1552 0.2225 155.7 3.460
. 4: 1 0.2 83.19 15.502 0.5817 0.8048 152.8 6.353
. 5: 1 0.3 75.88 21.354 1.2272 1.6382 149.6 8.752
. 6: 1 0.4 69.21 26.156 2.0463 2.6356 146.6 10.719
. 7: 1 0.5 63.13 30.046 3.0001 3.7278 144.1 12.314
. 8: 1 0.6 57.58 33.146 4.0553 4.8607 142.2 13.584
We need to know this to work the example; more details in a bit
Background
File name:
What would you like to “fix” in this plot?
0. [1] 0 6 12
model %>% intervention %>% simulate %>%
. # A tibble: 362 x 9
. ID time GUT CENT PERIPH CE EFFECT CP arm
. <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
. 1 1 0 0 0 0 0 157 0 1
. 2 1 0 100 0 0 0 157 0 1
. 3 1 0.1 91.2 8.44 0.155 0.222 156. 3.46 1
. 4 1 0.2 83.2 15.5 0.582 0.805 153. 6.35 1
. 5 1 0.3 75.9 21.4 1.23 1.64 150. 8.75 1
. 6 1 0.4 69.2 26.2 2.05 2.64 147. 10.7 1
. 7 1 0.5 63.1 30.0 3.00 3.73 144. 12.3 1
. 8 1 0.6 57.6 33.1 4.06 4.86 142. 13.6 1
. 9 1 0.7 52.5 35.6 5.18 5.99 141. 14.6 1
. 10 1 0.8 47.9 37.4 6.36 7.09 139. 15.3 1
. # … with 352 more rows
If we didn’t modify the output
. Model: effect
. Dim: 361 x 8
. Time: 0 to 36
. ID: 1
. ID time GUT CENT PERIPH CE EFFECT CP
. 1: 1 0.0 0 0 0 0 157 0
. 2: 1 0.1 0 0 0 0 157 0
. 3: 1 0.2 0 0 0 0 157 0
. 4: 1 0.3 0 0 0 0 157 0
. 5: 1 0.4 0 0 0 0 157 0
. 6: 1 0.5 0 0 0 0 157 0
. 7: 1 0.6 0 0 0 0 157 0
. 8: 1 0.7 0 0 0 0 157 0
mod %>% mrgsim() %>% mutate()mod %>% mrgsim() %>% filter()mod %>% mrgsim() %>% group_by()mod %>% mrgsim() %>% as_tibble()Demo: try 0.5 g Q 12 h for 1 day in normal and impaired renal function
File name:
On the fly
Persistent
mrgsim will call update for you (on the fly)
start, end, delta, addatol, rtolhmax, maxsteps, mxhnil, ixpr$OMEGA, $SIGMAtscale (rescale the output time)digits%>% %>% %>% model %>% intervention %>% %>% simulate %>% take-a-look
This is our second simulation workflow.
mod %>% ev(amt = 100, ii = 24, addl = 3) %>%idata_set( pop ) %>% mrgsim() %>% plot()
idata_set takes in individual-level data. ID CL VC KA KOUT IC50 FOO
. 1 1 1.05 47.8 0.8390 2.45 1.28 4
. 2 2 0.73 30.1 0.0684 2.51 1.84 6
. 3 3 2.82 23.8 0.1180 3.88 2.48 5
. [1] 10
. ID CL VC KA KOUT IC50 FOO
. 1 1 1.05 47.8 0.8390 2.45 1.28 4
. 2 2 0.73 30.1 0.0684 2.51 1.84 6
. 3 3 2.82 23.8 0.1180 3.88 2.48 5
.
. Model parameters (N=3):
. name value . name value
. CL 1 | V 20
. KA 1 | . .
idata?Batches of simulations or sensitivity analyses
. ID CL
. 1 1 0.50
. 2 2 0.75
. 3 3 1.00
. 4 4 1.25
. 5 5 1.50
+ mod %>%
idata_set(idata) %>%
ev(amt = 100, ii = 24, addl = 2) %>%
mrgsim(end = 120) %>%
plot(CP ~ ., logy=TRUE)Description: sensitivity analysis in a PBPK model
File name:
data_set is the dosing equivalent to idata_set. ID time amt ii addl cmt evid
. 1 1 0 100 24 3 1 1
. 2 2 0 300 24 3 1 1
. 3 3 0 1000 24 3 1 1
data_set configuration. ID time amt ii addl cmt evid CL
. 1 1 0 100 24 3 1 1 0.5
. 2 2 0 300 24 3 1 1 0.5
. 3 3 0 1000 24 3 1 1 0.5
. 4 4 0 100 24 3 1 1 1.0
. 5 5 0 300 24 3 1 1 1.0
. 6 6 0 1000 24 3 1 1 1.0
data <- mutate(data, dose = amt)
mod %>%
data_set(data) %>%
mrgsim(carry_out = "dose", end = 120) %>%
plot(log(CP)~time|factor(dose), group = ID, scales = "same")Copy from the input data to the output data
. Model: pk1
. Dim: 2886 x 4
. Time: 0 to 192
. ID: 3
. ID time dose CP
. 1: 2 0.0 300 0.000
. 2: 2 0.0 300 0.000
. 3: 2 0.2 300 2.712
. 4: 2 0.4 300 4.919
. 5: 2 0.6 300 6.712
. 6: 2 0.8 300 8.167
. 7: 2 1.0 300 9.345
. 8: 2 1.2 300 10.296
. ID time amt ii addl cmt evid CL dose
. 1 1 0 100 24 3 1 1 0.5 100
. 2 2 0 300 24 3 1 1 0.5 300
. 3 3 0 1000 24 3 1 1 0.5 1000
. Events:
. time amt ii addl cmt evid
. 1 0 100 24 3 1 1
mod + ev = ?mod + idata_set + ev = ?mod + data_set = ?. ID time amt ii addl cmt evid CL dose
. 1 1 0 100 24 3 1 1 0.5 100
. 2 2 0 300 24 3 1 1 0.5 300
. 3 3 0 1000 24 3 1 1 0.5 1000
. 4 4 0 100 24 3 1 1 1.0 100
. 5 5 0 300 24 3 1 1 1.0 300
. 6 6 0 1000 24 3 1 1 1.0 1000
Recall that data sets are just plain old data.frames in R. Feel free to code these up in any way that is convenient for you.
In our experience, the following helper functions cover many (not every) common needs for building the data sets.
expand.evas_data_setev_repev_assignexpand.evexpand.grid … gives “all combinations”. ID time amt cmt evid
. 1 1 0 100 1 1
. 2 2 0 100 1 1
. 3 3 0 200 1 1
. 4 4 0 200 1 1
as_data_setdata <- as_data_set(
ev(amt = 100, ii = 12, addl = 19, ID = 1:2),
ev(amt = 200, ii = 24, addl = 9, ID = 1:3),
ev(amt = 150, ii = 24, addl = 9, ID = 1:4)
). ID time cmt evid amt ii addl
. 1 1 0 1 1 100 12 19
. 2 2 0 1 1 100 12 19
. 3 3 0 1 1 200 24 9
. 4 4 0 1 1 200 24 9
. 5 5 0 1 1 200 24 9
. 6 6 0 1 1 150 24 9
. 7 7 0 1 1 150 24 9
. 8 8 0 1 1 150 24 9
. 9 9 0 1 1 150 24 9
ev_rep. time amt ii addl cmt evid ID
. 1 0 100 12 14 1 1 1
. 2 0 100 12 14 1 1 2
. 3 0 100 12 14 1 1 3
. 4 0 100 12 14 1 1 4
. 5 0 100 12 14 1 1 5
ev_rep. Events:
. time amt ii addl cmt evid
. 1 0 100 0 0 1 1
. 2 36 50 24 4 1 1
. time amt ii addl cmt evid ID
. 1 0 100 0 0 1 1 1
. 2 36 50 24 4 1 1 1
. 3 0 100 0 0 1 1 2
. 4 36 50 24 4 1 1 2
ev_assign. ID GROUP
. 1 1 1
. 2 2 2
. 3 3 1
. 4 4 2
. 5 5 1
. 6 6 2
e1 <- ev(amt = 100, ii = 24, addl = 9)
e2 <- mutate(e1, amt = 200)
data <- ev_assign(
l = list(e1,e2),
idata = pop,
evgroup = "GROUP",
join = TRUE
)
head(data). time amt ii addl cmt evid ID GROUP
. 1 0 100 24 9 1 1 1 1
. 2 0 200 24 9 1 1 2 2
. 3 0 100 24 9 1 1 3 1
. 4 0 200 24 9 1 1 4 2
. 5 0 100 24 9 1 1 5 1
. 6 0 200 24 9 1 1 6 2
Description: simulate GCSF PK/PD profiles for weight-based dosing
File name:
What’s going to happen?
What’s going to happen?
Combine two events
. Events:
. time amt cmt evid ii addl
. 1 0 200 1 1 0 0
. 2 24 100 1 1 24 2
What’s going to happen?
What’s going to happen?
Put two events in a sequence
. Events:
. time amt ii addl cmt evid
. 1 0 200 12 2 1 1
. 2 36 100 24 4 1 1
What is going to happen?
What is going to happen?
Wait before starting the next part of the regimen
. Events:
. time amt ii addl cmt evid
. 1 0 200 0 0 1 1
. 2 36 100 24 2 1 1
File name:
Section name:
%>% %>% %>%
mread or mread_cachemodlib("<model-name>")param(mod)init(mod)see(mod)ev(...): amt, cmt, time, ii, addl, rateDescription: Population PK of azithromycin
File name:
Description: simulate PK profiles using empirical Bayes estimates from a meropenem model fit in NONMEM
File name:
$PARAM [R] declare and initialize model parametersname = value format, or <newline>name except for words in mrgsolve:::reserved()R parsername anywhere in the modelR but read-only in the model specification file
$FIXED if appropriate$PARAM with the names in your input data sets$CMT [txt] declare model compartmentsname and number of compartments in the modelmrgsolve:::reserved()$INITExample $CMT
Example $INIT
$ODE [C++] write differential equationsdxdt_CMT = ratein - rateout;CMT for compartment amounts, parameters, or other variablesnew variables in your modelC++, which is a language$MAIN, $ODE, $TABLE, $GLOBALTo get a numeric variable, use double
Other types you might use
If in doubt, use double; it’s what you want most of the time.
$MAIN [C++] covariate model & initialsFor example
In this example
TVCL, KIN, KOUT, WT, SEX were declared in $PARAMRESPC++ expressions and functionsif(a == 2) b = 2;
if(b <= 2) {
c=3;
} else {
c=4;
}
double d = pow(base,exponent);
double e = exp(3);
double f = fabs(-4);
double g = sqrt(5);
double h = log(6);
double i = log10(7);
double j = floor(4.2);
double k = ceil(4.2);Lots of help on the web http://en.cppreference.com/w/cpp/numeric/math/tgamma
#define)POPULATION simulation$OMEGA and $SIGMA to define covariance matrices for IIV and RUV, respectively
ETA(n) and EPS(m) for realized random effects drawn from $OMEGA and $SIGMA respectively
ETAs insteadmrgsolve recognizes ID as a subject indicator for simulating a new ETA
mrgsolve allows you to write an error model as a function of EPS(m) and any other calculated value in your model
$OMEGA and $SIGMA [txt]$OMEGA and $SIGMA are allowed and each may be namedmrgsolve will enforce later on$OMEGA and $SIGMA@ delimited & must be on different line from the matrix data@block)$OMEGA and $SIGMA$OMEGA and $SIGMA are allowed; each may be namedUsers are encouraged to add labels
@ macros$OMEGA and $SIGMA@ should appear at the start of the line and the entire line is reserved for options@cor @name PKTwo forms:
@block means block=TRUE@labels ECL EVC means labels=c("ECL", "EVC")$PKMODELncmt to 1 or 2depot to TRUE for extravascular dosing compartmenttrans to change the names of required parameters$CMT to declare 1 to 3 compartments as appropriate$PKMODELTo change the bioavability of doses administered to a compartment, set F_CMT in $MAIN
To add a lagtime for doses administered to a compartment, set ALAG_CMT in $MAIN
D_CMT in $MAIN
rate = -2 in your data set or event objectR_CMT in $MAIN
rate = -1 in your data set or event objectFirst time to read
Next time to read
soloc)It can be helpful to set
and
We
. ID time GUT CENT CP
. 1 1 0 0 0 0
. 2 1 1 0 0 0
. ID time CP
. 1 1 0 0
. 2 1 1 0
. ID time GUT CENT CP
. 1 1 0 0 0 0
. 2 1 0 1 0 0
. ID time GUT CENT CP
. 1 1 0 0.0000000 0.0000000 0.00000000
. 2 1 1 0.3011942 0.6782976 0.03391488